This script analyzes and plots data for Symbiontic Integration 2021 respirometry data for thermal performance curves.
Set up workspace, set options, and load required packages.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse')
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car')
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales')
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2')
if ("nls.multstart" %in% rownames(installed.packages()) == 'FALSE') install.packages('nls.multstart')
if ("broom" %in% rownames(installed.packages()) == 'FALSE') install.packages('broom')
if ("rTPC" %in% rownames(installed.packages()) == 'FALSE') remotes::install_github("padpadpadpad/rTPC")
#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('scales')
library('nls.multstart')
library('broom')
library('rTPC')
library('ggstatsplot')
Load data from LoLinR.
PRdata<-read.csv("Pacu2021/Output/Respiration/TPC/oxygen_Resp_rates.csv") #load data
Format data columns.
#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]
#format columns
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Temperature<-as.numeric(as.character(PRdata$Temperature))
PRdata$Treatment<-as.factor(PRdata$Treatment)
PRdata$Plate<-as.factor(PRdata$Plate)
PRdata$Run<-as.factor(PRdata$Run)
PRdata <- PRdata %>%
group_by(Temperature,Treatment)%>%
mutate(abs=abs(R.nmol.org.min))%>%
summarise(rate=abs) %>%
mutate(temp=Temperature)
PRdata$group <- paste0(PRdata$Temperature,"_", PRdata$Treatment)
View and remove outliers.
#PRdata <- PRdata %>% filter(R.nmol.org.min > -0.15) %>% filter(R.nmol.org.min < 0)
#identify outliers by Tremperature and Treatment groups
outlier.plot <- ggbetweenstats(PRdata,group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).
outlier.plot
<<<<<<< HEAD
ggsave("Pacu2021/Figures/Respiration/TPC/outliersbygroup.pdf", outlier.plot, dpi=300, w=16, h=8, units="in")
q <- c(0.25, 0.75)
Quants <- PRdata %>%
group_by(Treatment, Temperature) %>%
summarize(quant25 = quantile(rate, probs = q[1]),
quant75 = quantile(rate, probs = q[2]),
IQRbyGroup=IQR(rate))
Quants$group <-paste0(Quants$Temperature,"_", Quants$Treatment)
Quants$upper <- Quants$quant75+1.5*Quants$IQRbyGroup # Upper Range
Quants$lower <- Quants$quant25-1.5*Quants$IQRbyGroup # Lower Range
#join outlier cutoffs with data
PRdata <- left_join(PRdata, Quants)
#remove outliers
PRdata <- PRdata %>%
filter(rate < upper) %>%
filter(rate > lower)
outlier.removed.plot <- ggbetweenstats(PRdata, group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).
outlier.removed.plot
<<<<<<< HEAD
ggsave("Pacu2021/Figures/Respiration/TPC/Removed_outliersbygroup.pdf", outlier.removed.plot, dpi=300, w=16, h=8, units="in")
Plot data as a dot plot
r_plot<-PRdata %>%
ggplot(., aes(x = Temperature, y = rate, colour=Treatment)) +
geom_point(aes(fill=Treatment, group=interaction(Temperature, Treatment)), pch = 21, size=2, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Treatment", values=c("blue","red"))+
scale_color_manual(name="Treatment", values=c("blue","red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot
<<<<<<< HEAD
ggsave("Pacu2021/Figures/Respiration/TPC/Respiration_Scatter.pdf", r_plot, dpi=300, w=8, h=8, units="in")
Filtering
PRdata <- PRdata %>%
filter(!Temperature==8) %>%
filter(!Temperature==18)%>%
filter(!Temperature==34)%>%
filter(!Temperature==36)%>%
filter(!Temperature==38)%>%
filter(!Temperature==40)
r_plot<-PRdata %>%
ggplot(., aes(x = Temperature, y = rate, colour=Treatment)) +
geom_point(aes(fill=Treatment, group=interaction(Temperature, Treatment)), pch = 21, size=2, alpha=0.5, position = position_jitterdodge(0.3)) +
#geom_smooth(method = "lm", se=FALSE, formula=y ~ poly(x, 3, raw=TRUE), size=3)+
xlab("Temperature") +
scale_fill_manual(name="Treatment", values=c("blue","red"))+
scale_color_manual(name="Treatment", values=c("blue","red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot
<<<<<<< HEAD
ggsave("Pacu2021/Figures/Respiration/TPC/Respiration_Scatter_Filtered.pdf", r_plot, dpi=300, w=8, h=8, units="in")
View and remove outliers.
#identify outliers by Tremperature and Treatment groups
outlier.plot <- ggbetweenstats(PRdata,group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).
outlier.plot
<<<<<<< HEAD
Quants <- PRdata %>%
group_by(Treatment, Temperature) %>%
summarize(quant25 = quantile(rate, probs = q[1]),
quant75 = quantile(rate, probs = q[2]),
IQRbyGroup=IQR(rate))
Quants$group <-paste0(Quants$Temperature,"_", Quants$Treatment)
Quants$upper <- Quants$quant75+1.5*Quants$IQRbyGroup # Upper Range
Quants$lower <- Quants$quant25-1.5*Quants$IQRbyGroup # Lower Range
#join outlier cutoffs with data
PRdata <- left_join(PRdata, Quants)
#remove outliers
PRdata <- PRdata %>%
filter(rate < upper) %>%
filter(rate > lower)
<<<<<<< HEAD
PRdata <- PRdata %>% filter(rate > 0.0015)
=======
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
PRdata <- PRdata %>% filter(rate < 0.14)
PRdata <- PRdata[-c(96,106),]
outlier.removed.plot <- ggbetweenstats(PRdata, group, rate, outlier.tagging = TRUE)
## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).
outlier.removed.plot
<<<<<<< HEAD
ggsave("Pacu2021/Figures/Respiration/TPC/Removed_outliersbygroup2.pdf", outlier.removed.plot, dpi=300, w=16, h=8, units="in")
#Padifeld et al rTPC and nls.multstart: A new pipeline to fit thermal performance curves in r ### https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13585
# choose model
get_model_names()
## [1] "beta_2012" "boatman_2017" "briere2_1999"
## [4] "delong_2017" "flinn_1991" "gaussian_1987"
## [7] "hinshelwood_1947" "joehnk_2008" "johnsonlewin_1946"
## [10] "kamykowski_1985" "lactin2_1995" "lrf_1991"
## [13] "modifiedgaussian_2006" "oneill_1972" "pawar_2018"
## [16] "quadratic_2008" "ratkowsky_1983" "rezende_2019"
## [19] "sharpeschoolfull_1981" "sharpeschoolhigh_1981" "sharpeschoollow_1981"
## [22] "spain_1982" "thomas_2012" "thomas_2017"
## [25] "weibull_1995"
#sharpeschoolhigh_1981
# get start vals
start_vals <- get_start_vals(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')
# get limits
low_lims <- get_lower_lims(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')
upper_lims <- get_upper_lims(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')
start_vals
<<<<<<< HEAD
## r_tref e eh th
## 0.0456096 0.7027915 2.1208907 30.9411765
=======
## r_tref e eh th
## 0.04434088 0.76808284 2.42820170 30.94117647
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
low_lims
## r_tref e eh th
## 0 0 0 1
upper_lims
## r_tref e eh th
## 0.1232836 10.0000000 20.0000000 32.0000000
# AMBIENT CURVE FIT
d.amb <- PRdata %>%
filter(Treatment=="Ambient")
# # get start vals
# amb.start_vals <- get_start_vals(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.start_vals
#
# # get limits
# amb.low_lims <- get_lower_lims(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.low_lims
# amb.upper_lims <- get_upper_lims(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.start_vals
# amb.low_lims
# amb.upper_lims
amb.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.amb,
iter = 500,
start_lower = start_vals - 1,
start_upper = start_vals + 1,
lower = low_lims,
upper = upper_lims,
supp_errors = 'Y')
amb.fit
## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
<<<<<<< HEAD
## 0.06369 0.77099 8.81372 31.66205
## residual sum-of-squares: 0.03744
##
## Number of iterations to convergence: 25
=======
## 0.06287 0.81631 10.05527 31.49888
## residual sum-of-squares: 0.03564
##
## Number of iterations to convergence: 21
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
## Achieved convergence tolerance: 1.49e-08
# predict new data
amb_new_data <- data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), 0.5))
amb.preds <- augment(amb.fit, newdata = amb_new_data)
amb.TCP.res <- calc_params(amb.fit) %>%
mutate_all(round, 2) # round
amb.TCP.res
<<<<<<< HEAD
## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.07 29.55 0.81 34.96 0.7 2.31 2.5 5.41 34.15
## breadth skewness
## 1 4.99 -1.61
=======
## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.07 29.58 2.34 34.35 0.77 2.9 2.71 4.77 32.01
## breadth skewness
## 1 4.58 -2.14
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
Resampling the original data with replacement
# refit model using nlsLM
amb.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.amb,
start = coef(amb.fit),
lower = low_lims,
upper = upper_lims,
weights = rep(1, times = nrow(d.amb)))
# bootstrap using case resampling
amb.boot1 <- Boot(amb.fit_nlsLM, method = 'case')
# look at the data
head(amb.boot1$t)
## r_tref e eh th
<<<<<<< HEAD
## [1,] 0.05622679 1.0115372 3.274674 32.00000
## [2,] 0.06621550 0.6028136 12.494659 31.84514
## [3,] 0.11625218 1.5155506 4.932409 29.19690
## [4,] 0.12328360 1.6623770 3.525770 28.08559
## [5,] 0.12328360 1.5686641 2.828922 27.31338
## [6,] 0.10076807 1.2078562 3.945774 29.13973
=======
## [1,] 0.12328360 1.8657450 4.027641 28.04144
## [2,] 0.09298738 2.4680678 7.536008 29.25540
## [3,] 0.06189850 0.8984202 20.000000 31.68621
## [4,] 0.12328360 2.0854617 5.215647 28.79464
## [5,] 0.07015725 1.0264364 20.000000 31.58002
## [6,] 0.07430899 1.1065485 9.812652 31.17470
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
# create predictions of each bootstrapped model
amb.boot1_preds <- amb.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
amb.boot1_conf_preds <- group_by(amb.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
amb.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
amb.CI.plot
<<<<<<< HEAD
# HIGH CURVE FIT
d.high <- PRdata %>%
filter(Treatment=="High")
# choose model
mod <- 'sharpschoolhigh_1981'
# # get start vals
# high.start_vals <- get_start_vals(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.start_vals
#
# # get limits
# high.low_lims <- get_lower_lims(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.low_lims
# high.upper_lims <- get_upper_lims(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.start_vals
# high.low_lims
# high.upper_lims
high.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.high,
iter = 500,
start_lower = start_vals - 1,
start_upper = start_vals + 1,
lower = low_lims,
upper = upper_lims,
supp_errors = 'Y')
high.fit
## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
## 0.1233 1.4038 2.4282 26.4288
## residual sum-of-squares: 0.03235
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08
high_new_data <- data.frame(temp = seq(min(d.high$temp), max(d.high$temp), 0.5))
high.preds <- augment(high.fit, newdata = high_new_data)
high.TCP.res <- calc_params(high.fit) %>%
mutate_all(round, 2) # round for easy viewing
high.TCP.res
## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.44 7.45 57.05 0.53 0.47 2.01 29.62 49.6
## breadth skewness
## 1 8.87 0.06
Resampling the original data with replacement
# refit model using nlsLM
high.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.high,
start = coef(high.fit),
lower = low_lims,
upper = upper_lims,
weights = rep(1, times = nrow(d.high)))
# bootstrap using case resampling
<<<<<<< HEAD
high.boot1 <- Boot(high.fit_nlsLM, method = 'case')
# look at the data
head(high.boot1$t)
## r_tref e eh th
## [1,] 0.1232836 1.0875773 1.3645045 24.82696
## [2,] 0.1232836 1.0663405 1.7193422 25.18904
## [3,] 0.1232836 1.5667784 3.4832088 27.14641
## [4,] 0.1232836 1.2804056 2.9616468 26.35147
## [5,] 0.1232836 1.5242106 3.1892471 27.27148
## [6,] 0.1232836 0.0798746 0.1515831 1.00000
=======
high.boot1 <- Boot(high.fit_nlsLM, method = 'case')
##
## Number of bootstraps was 997 out of 999 attempted
# look at the data
head(high.boot1$t)
## r_tref e eh th
## [1,] 0.09466776 1.1017541 3.245382 28.74252
## [2,] 0.12328360 1.5033158 2.816571 26.55534
## [3,] 0.12328360 1.3592890 2.435695 26.91997
## [4,] 0.12328360 1.6925379 2.945170 26.25989
## [5,] 0.07216997 0.8133763 4.153968 29.78028
## [6,] 0.04920720 0.3859032 4.080503 32.00000
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
# create predictions of each bootstrapped model
high.boot1_preds <- high.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.high$temp), max(d.high$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
high.boot1_conf_preds <- group_by(high.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
high.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), high.preds, col = 'red') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
high.CI.plot
<<<<<<< HEAD
# plot data and model fit
TPC.plot <- ggplot() +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue', size=2) +
geom_line(aes(temp, .fitted), high.preds, col = 'red', size=2) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Metabolic rate')
TPC.plot
<<<<<<< HEAD
ggsave("Pacu2021/Figures/Respiration/TPC/TPC.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")
=======
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_SharpSchool.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
#confidence intervals of TPC parameters
broom::tidy(amb.fit_nlsLM)
## # A tibble: 4 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
<<<<<<< HEAD
## 1 r_tref 0.0637 0.0110 5.78 5.07e- 7
## 2 e 0.771 0.412 1.87 6.76e- 2
## 3 eh 8.81 7.46 1.18 2.43e- 1
## 4 th 31.7 0.850 37.2 1.34e-37
=======
## 1 r_tref 0.0629 0.00968 6.50 3.67e- 8
## 2 e 0.816 0.385 2.12 3.89e- 2
## 3 eh 10.1 7.32 1.37 1.75e- 1
## 4 th 31.5 0.755 41.7 1.56e-40
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
broom::tidy(high.fit_nlsLM)
## # A tibble: 4 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 r_tref 0.123 0.497 0.248 0.805
## 2 e 1.40 3.76 0.373 0.710
## 3 eh 2.43 1.14 2.12 0.0384
## 4 th 26.4 20.7 1.28 0.208
#AMBIENT
amb.extra_params <- calc_params(amb.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
amb.ci_extra_params <- Boot(amb.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(amb.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
##
<<<<<<< HEAD
## Number of bootstraps was 83 out of 200 attempted
=======
## Number of bootstraps was 68 out of 200 attempted
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
amb.ci_extra_params <- left_join(amb.ci_extra_params, amb.extra_params)
amb.ci_extra_params$Treatment <- "Ambient"
# Joining, by = "param"
#HIGH
high.extra_params <- calc_params(high.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
high.ci_extra_params <- Boot(high.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(high.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
##
<<<<<<< HEAD
## Number of bootstraps was 139 out of 200 attempted
=======
## Number of bootstraps was 151 out of 200 attempted
high.ci_extra_params <- left_join(high.ci_extra_params, high.extra_params)
high.ci_extra_params$Treatment <- "High"
#Join Ambient and High estimates and CIs
All_params <- rbind(amb.ci_extra_params, high.ci_extra_params)
All_params <- All_params %>%
mutate_if(is.numeric, round, 2)
#Plot
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
geom_point(size = 2) +
scale_color_manual(name="Treatment", values=c("blue","red"))+
geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
theme_bw() +
facet_wrap(~param, scales = 'free_y') +
scale_x_discrete('')
estimate.plots
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates_SharpSchool.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")
#Gaussian
#AMBIENT
amb.gaus.fit <- nls_multstart(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.amb,
iter = c(4,4,4),
start_lower = get_start_vals(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987') - 1,
start_upper = get_start_vals(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987') + 1,
lower = get_lower_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
supp_errors = 'Y',
convergence_count = FALSE)
amb.gaus.fit
## Nonlinear regression model
## model: rate ~ gaussian_1987(temp = temp, rmax, topt, a)
## data: data
## rmax topt a
## 0.06347 27.83799 4.38800
## residual sum-of-squares: 0.03694
##
## Number of iterations to convergence: 11
## Achieved convergence tolerance: 1.49e-08
# predict new data
amb_new_data <- data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), 0.5))
amb.preds <- augment(amb.gaus.fit, newdata = amb_new_data)
amb.TCP.res <- calc_params(amb.gaus.fit) %>%
mutate_all(round, 2) # round for easy viewing
amb.TCP.res
## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.06 27.84 17.1 38.58 1.37 0.73 NA 10.74 21.48
## breadth skewness
## 1 5.87 0.65
amb.gaus.fit_nlsLM <- minpack.lm::nlsLM(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.amb,
start = coef(amb.gaus.fit),
lower = get_lower_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
weights = rep(1, times = nrow(d.amb)))
# bootstrap using case resampling
amb.boot1 <- Boot(amb.gaus.fit_nlsLM, method = 'case')
# look at the data
head(amb.boot1$t)
## rmax topt a
## [1,] 0.06197087 27.13551 4.854462
## [2,] 0.05900555 27.60181 4.673603
## [3,] 0.05973728 28.41771 4.843638
## [4,] 0.06388289 27.78184 4.564628
## [5,] 0.05356135 27.06235 5.075532
## [6,] 0.07453052 28.19593 3.481052
# create predictions of each bootstrapped model
amb.boot1_preds <- amb.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = gaussian_1987(temp = temp, rmax, topt, a))
# calculate bootstrapped confidence intervals
amb.boot1_conf_preds <- group_by(amb.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
amb.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
amb.CI.plot
# HIGH
high.gaus.fit <- nls_multstart(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.high,
iter = c(4,4,4),
start_lower = get_start_vals(d.high$temp, d.high$rate, model_name = 'gaussian_1987') - 1,
start_upper = get_start_vals(d.high$temp, d.high$rate, model_name = 'gaussian_1987') + 1,
lower = get_lower_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
supp_errors = 'Y',
convergence_count = FALSE)
high.gaus.fit
## Nonlinear regression model
## model: rate ~ gaussian_1987(temp = temp, rmax, topt, a)
## data: data
## rmax topt a
## 0.04746 27.07125 5.77867
## residual sum-of-squares: 0.0274
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08
# predict new data
high_new_data <- data.frame(temp = seq(min(d.high$temp), max(d.high$temp), 0.5))
high.preds <- augment(high.gaus.fit, newdata = high_new_data)
high.TCP.res <- calc_params(high.gaus.fit) %>%
mutate_all(round, 2) # round for easy viewing
high.TCP.res
## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.07 12.93 41.22 0.53 0.79 2.01 14.15 28.29
## breadth skewness
## 1 7.72 -0.26
high.gaus.fit_nlsLM <- minpack.lm::nlsLM(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.high,
start = coef(high.gaus.fit),
lower = get_lower_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
weights = rep(1, times = nrow(d.high)))
# bootstrap using case resampling
high.boot1 <- Boot(high.gaus.fit_nlsLM, method = 'case')
# look at the data
head(high.boot1$t)
## rmax topt a
## [1,] 0.05927776 26.86467 4.041476
## [2,] 0.05146475 26.79550 4.483799
## [3,] 0.04765019 30.66779 8.565635
## [4,] 0.04808687 28.83889 8.571844
## [5,] 0.04426714 25.74916 8.505839
## [6,] 0.04689024 28.19428 6.521274
# create predictions of each bootstrapped model
high.boot1_preds <- high.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.high$temp), max(d.high$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = gaussian_1987(temp = temp, rmax, topt, a))
# calculate bootstrapped confidence intervals
high.boot1_conf_preds <- group_by(high.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
high.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), high.preds, col = 'red') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
high.CI.plot
# plot data and model fit
TPC.plot <- ggplot() +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue', size=2) +
geom_line(aes(temp, .fitted), high.preds, col = 'red', size=2) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Metabolic rate')
TPC.plot
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_Gaussian.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")
amb.TCP.res
## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.06 27.84 17.1 38.58 1.37 0.73 NA 10.74 21.48
## breadth skewness
## 1 5.87 0.65
high.TCP.res
## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.07 12.93 41.22 0.53 0.79 2.01 14.15 28.29
## breadth skewness
## 1 7.72 -0.26
#confidence intervals of Gaussian TPC parameters
broom::tidy(amb.gaus.fit_nlsLM)
## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 rmax 0.0635 0.00604 10.5 2.38e-14
## 2 topt 27.8 0.570 48.9 1.67e-44
## 3 a 4.39 0.778 5.64 7.48e- 7
broom::tidy(high.gaus.fit_nlsLM)
## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 rmax 0.0475 0.00492 9.65 2.91e-13
## 2 topt 27.1 0.847 32.0 2.56e-36
## 3 a 5.78 1.60 3.61 6.71e- 4
#AMBIENT
amb.extra_params <- calc_params(amb.gaus.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
amb.ci_extra_params <- Boot(amb.gaus.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(amb.gaus.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
##
## Number of bootstraps was 109 out of 200 attempted
amb.ci_extra_params <- left_join(amb.ci_extra_params, amb.extra_params)
amb.ci_extra_params$Treatment <- "Ambient"
# Joining, by = "param"
#HIGH
high.extra_params <- calc_params(high.gaus.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
high.ci_extra_params <- Boot(high.gaus.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(high.gaus.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')
##
## Number of bootstraps was 159 out of 200 attempted
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
high.ci_extra_params <- left_join(high.ci_extra_params, high.extra_params)
high.ci_extra_params$Treatment <- "High"
#Join Ambient and High estimates and CIs
All_params <- rbind(amb.ci_extra_params, high.ci_extra_params)
All_params <- All_params %>%
mutate_if(is.numeric, round, 2)
#Plot
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
geom_point(size = 2) +
scale_color_manual(name="Treatment", values=c("blue","red"))+
geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
theme_bw() +
facet_wrap(~param, scales = 'free_y') +
scale_x_discrete('')
estimate.plots
<<<<<<< HEAD
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")
=======
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates_Gaussian.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695